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ABSTRACT 

Galaxies and the dark matter halos that host them are not spherically symmetric, yet spherical 
symmetry is a helpful simplifying approximation for idealised calculations and analysis of 
observational data. The assumption leads to an exact conservation of angular momentum for 
every particle, making the dynamics unrealistic. But how much does that inaccuracy matter in 
practice for analyses of stellar distribution functions, collisionless relaxation, or dark matter 
core-creation? 

We provide a general answer to this question for a wide class of aspherical systems; 
specifically, we consider distribution functions that are “maximally stable”, i.e. that do not 
evolve at first order when external potentials (which arise from baryons, large scale tidal 
fields or infalling substructure) are applied. We show that a spherically-symmetric analysis of 
such systems gives rise to the false conclusion that the density of particles in phase space is 
ergodic (a function of energy alone). 

Using this idea we are able to demonstrate that: (a) observational analyses that falsely 
assume spherical symmetry are made more accurate by imposing a strong prior preference 
for near-isotropic velocity dispersions in the centre of spheroids; (b) numerical simulations 
that use an idealised spherically-symmetric setup can yield misleading results and should be 
avoided where possible; and (c) triaxial dark matter halos (formed in collisionless cosmologi¬ 
cal simulations) nearly attain our maximally-stable limit, but their evolution freezes out before 
reaching it. 


1 INTRODUCTION 


Spherical symmetry is a foundational assumption of many dynami¬ 
cal analyses. The primary motivation is simplicity, since few astro¬ 
nomical objects are actually spherical. For example, observations 
and simulations both suggest that gravitational potential wells gen¬ 


erated by dark matter halos are typically triaxial (e.g. Dubinski & 


Carlberg 1991 

Cole & Lacey 1996 

|Jing & Suto 2002; Kasun & 

Evrard 2005; 

layashi et al. 2007 

Schneider et al. 2012 Loeb- 


averaged densities and velocities (e.g. Dubinski & Carlberg 1991 
|Navarro et al.||1996{ |Taylor & Navarro||2001[ Stadel et al. 2009| 

at best tells only part of the story. At worst, it could be severely 
misleading. 


The question of whether baryonic processes can convert dark 
matter cusps into cores (Pontzen & Governato |2014[ l provides one 
motivation for a detailed study of the relationship between spherical 


and near-spherical dynamics. To explain why, we need to look for¬ 
ward to some of our results. Later in this paper, we cut a dark matter 
halo out of a cosmological simulation, then match it to an exactly 
spherical halo with an identical density and velocity anisotropy pro¬ 
file. This gives us two easy-to-compare equilibrium structures - the 
first triaxial, the second spherical - to perform a dynamical com¬ 
parison. We expose each to the same time-varying gravitational po¬ 
tential, mimicking the effects of stellar feedback (there are no ac¬ 
tual baryons in these runs). After 1 Gyr, the triaxial halo’s averaged 
density profile flattens into a convincing dark matter core, but the 
spherical halo maintains its cusp (see Figure 1). 

This example, which is fully explored in Section [3~5| illus¬ 
trates how it is dangerous to use spherically-symmetric simulations 
to infer anything about dynamics — even spherically-averaged dy¬ 
namics — in the real universe. A spherical system does not evolve 
in the same way as the spherical averages of a triaxial system. 

Ignoring asphericity can also lead to observational biases (e.g. 
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2 A. Pontzen et al 


|Hayashi & Navarro|2006[|Corless & King|2007j i. From a dynamical 
standpoint, the nature of orbits in triaxial potentials is fundamen¬ 
tally different from those in spherical potentials: although the to¬ 
tal angular momentum of any self-gravitating system must always 
be conserved, it is only in the spherical case that this conserva¬ 
tion holds for individual particles. Conversely, a large fraction of 
dark matter particles near the centre of cosmological halos will be 
on box orbits which do not conserve their individual angular mo¬ 


menta i deZeeuw&Memtt|1983^ Merritt & Valluri|1996[ |Holley-| 


|Bockelmann et al.||2001 Adams et al.||2007|. One practical 


sequence is that asphericity may be responsible for filling the loss 
cones of supermassive black holes at the centre of the correspond¬ 
ing galaxies jMerritt & Poon|2004| >. 

Finally, it is known that asphericity plays a fundamental role 
in setting the equilibrium density profile during gravitational cold 


as the radial orbit instability or ROI 

Henon 1973; Saffl 1991 Huss 

et al. 1999; 

MacMillan et al. 2006 

Bellovary et al. 20081 Barnes 

et al. 2009 [ Marechal & Perez 2009 

; a related effect was discussed 

by Adams et al. (2007). The name 

arises because particles on ra- 


dial orbits are perturbed onto more circular trajectories. At the same 
time, the density distribution becomes triaxial. Even in the case of a 
uniform spherical collapse, this symmetry-breaking process is still 
triggered, presumably by numerical noise; the tangential compo¬ 
nent of forces must be unphysically suppressed for the system to 
remain spherical ( |Huss et al.|1999[|MacMillan et al.|2006| . 

Despite all this, assuming spherical symmetry is very tempting 
because it makes life so much easier. Defining spherically-averaged 
quantities is a well-defined and sensible procedure even if we actu¬ 
ally have the full distribution function in hand (as in simulations): 
departures from spherical symmetry are sufficiently small that dif¬ 
ferent averaging procedures lead to consistent results jSaha & Read| 
|2009[ >. Additionally, when an aspherical halo is in equilibrium, we 
have shown numerically that a “sphericalised” version of it is also 
in equilibrium (see Appendix B of[Pontzen & Governato 2013). 
This is helpful because it allows one to make a meaningful analysis 
in spherical coordinates, even when the system is aspherical. But it 
breaks down when out-of-equilibrium processes are included, as in 
the stellar-feedback-driven core-creation example above. 

The present paper formalizes the idea of spherical analysis 
performed on aspherical systems and follows through the conse¬ 
quences. We will study equilibrium distribution functions in nearly, 
but not exactly, spherically-symmetric potentials, and focus on 
maximally-stable systems (which we define as being stable against 
all possible external linear perturbations). We will find that, in 
spherical coordinates, such systems appear to be “ergodic'’ (mean¬ 
ing that their distribution functions depend on energy alone) be¬ 
cause the individual particles move randomly in angular momen¬ 
tum while maintaining a near-constant energy. It is important to 
emphasise that this describes the appearance of the system when 
analysed in a spherical coordinate system and the true system need 
not be chaotic for the result to hold, provided any isolating integrals 
are not closely related to angular momentum. 

The formal statement of this idea is derived in Section [2] A 
brief overview of the required background is given in Appendix [A] 
and a second-order derivation of the evolution is given in Ap¬ 
pendix [b] Section [3] outlines the practical consequences, starting 
by recasting and extending the phenomenology of the radial or¬ 
bit instability. We describe an immediate implication for observa¬ 
tional studies of aspherical systems, namely a new way to break the 
anisotropy degeneracy. Finally we return to the motivating problem 
above and explain why triaxial systems can undergo cusp-core tran- 



Figure 1 . One motivation for studying the relationship between spherical 
and aspherical dynamics is that the conversion of a dark matter cusp into a 
core by baryonic processes is qualitatively different in the two cases. Here 
we show the inner log density slope from a numerical experiment on two ha¬ 
los. One is spherical (dashed line) and the other triaxial (solid line) but their 
spherically-averaged properties are initially identical. An external potential 
has been added at the centre during the times indicated by the grey bands, 
with the fluctuations mimicking stellar feedback. The triaxial halo develops 
a clear core, whereas the spherical halo almost maintains its central density 
cusp. A complete description and analysis is given in Section|33| 


sitions more easily than spherical systems. Section[4]concludes and 
points to open questions and future work. 


2 ASPHERICAL DYNAMICS IN SPHERICAL 
COORDINATES 

In this section we consider an aspherical system which is maxi¬ 
mally stable against external linear perturbations. We assume that 
an observer of this system analyses it assuming spherical symmetry. 
We will show that this observer (falsely) concludes that the system 
is ergodic, i.e. that the density of particles in phase space is a func¬ 
tion of energy alone. The derivation requires the use of action-angle 
coordinates; a crash course is provided in Appendix [A] 


2.1 Single particles 

Given any near-spherical system, the Hamiltonian in the spherical 
action-angle variables is 

H(J,@) = Ho(J) + SH(J,&), (1) 

where Hq is the sum of kinetic and potential energies in the spher¬ 
ical background, J = ( J r ,j,jz ) is the vector of spherical actions 
(see Appendix [A|, © is the vector of spherical angles and 8H con¬ 
tains the perturbation (which includes the aspherical correction to 
the potential). The orbit of a particle in exact spherical symmetry, 
SH = 0, is described by Hamilton's equations: 

= flo(Jo). (2) 

J=Jo 

which defines the constant background orbital frequencies Qo(Jo)- 
The expressions Jq, ©0 and SIq(Jq) will be used throughout to re¬ 
fer to the background (SH = 0) solution. This algebraically simple 


Jn = - 


dHo 

d& 


= 0 ; 


©n = 


dH 0 

dJ 
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form of the equations of motion is the reason for using action-angle 
variables, since it immediately integrates to 

Jo(t) = Jo = constant; & 0 (t) = 0 o (O) + Sl o (Jo)t, (3) 


where Jo and ©y(0) specify the initial action and angle coordinates 
of the orbit. 

We now consider the effect of the aspherical correction to the 
potential encoded in SH, using standard Hamiltonian perturbation 
theory (e.g. Lichtenberg & Lieberman 1992, Bi nney & Tremaine| 
2008). First, taking advantage of the angle coordinates 0 being 
periodic in 2k, 8H is expressed as 

8H(J,G) = Y,8 H n(J)e in& . ( 4 ) 

n 

This equation states that, at any fixed J, one can expand the peri¬ 
odic 0 dependence in a Fourier series without loss of generality. 

We are interested in the evolution of J at first order in the 
perturbation. Hamilton’s relevant equation now reads: 

J = -^ = -'£in8H n (J)e in& . (5) 

d& n 

Because 8H is small, the result to first-order accuracy is given by 
substituting the zero-order solution 0 into equation 0 and inte¬ 
grating to give 

j(t) = Jo - L ™ SH rl ( Jo)e m & ri) +.... (6) 

™ n “o 

Consequently as n ■ S~2q —> 0, the linear-order correction to the orbit 
of a particle can become large even if the aspherical correction to 
the potential (SH) is small, an effect known as resonance (Binney 
|& Trem aine 20081. Consider now the evolution of the background 
energy along the perturbed trajectory, given by 

r) M 

Hq(J (f)) =* Ho(Jo) + -jj ■ (, J(t) - Jo) + • • •, (7) 

where we have Taylor-expanded to first order around Jq. Substitut¬ 
ing equation [6] for J(t), there is a cancellation between numerator 
and denominator: 


H 0 ( J(0) - Ho( Jo) ~ ^SH n (J 0 )e in & ^ + ■ 


( 8 ) 


and so the fractional variation in Ho remains small, even if J 
changes significantly over time. In other words, according to lin¬ 
ear perturbation theory, particles migrate large distances in J along 
surfaces of approximately constant background energy Hq. One can 
verify this constrained-migration prediction in numerical simula¬ 
tions of dark matter halos - an explicit demonstration is given in 
Appendix |A3| This gives us some intuition for the result to come: a 
population of particles will seem to “randomise” their actions (in¬ 
cluding angular momentum), but not their energy distribution. 

The extent of the migration will depend on the nature of the 
potential in which a particle orbits. To quantify this requires going 
beyond linear perturbation theory and is the subject of “KAM the¬ 
ory” after Kolmogorov ( 1 954}, |Arnold 1 1963 i and|Moser i) 1 962|i; 


see e.g. Binney & Tremaine (2008 1 ; Goldstein et al.j i 2002 ; Licht-| 
|enberg & Lt eberman (1992 ) for introductions. Colloquially the re¬ 
sult is that for any given small perturbation the migration of typical 
orbits is also small. Arnold diffusion offers the most famous route 
to more significant diffusion through action space (see[Lichtenberg| 
|& Lieb erman 1992J; but in our case, there is a more immediate rea¬ 
son why the KAM result does not in fact hold. Specifically, KAM 
theory relies on the frequencies fi(J) being non-degenerate - i.e. 
that any change in the action leads also to a change in the frequen¬ 
cies, thus shutting off resonant migration. In smooth potentials, O. 


is almost a function of energy alone (see Appendix |A3| l and so the 
migration can be substantial. 

Overall we informally expect particles to redistribute them¬ 
selves randomly within the action shell of fixed background energy 
until they are evenly spread, implying a distribution function that 
appears ergodic in a spherical analysis. This does not imply the 
orbits are chaotic in the traditional sense; it is only because we 
are analysing an aspherical system in spherical coordinates that the 
phenomenon arises. With this in mind, we now turn to a more for¬ 
mal demonstration of the result. 


2.2 The distribution function 


So far we have discussed how a single particle orbiting in a mildly 
aspherical potential does not conserve its spherical actions (e.g. an¬ 
gular momentum). We informally suggested that a population will 
appear to ‘randomise’ the spherical actions at fixed energy. We now 
show more formally that a distribution function of particles subject 
to aspherical perturbations will be most stable when it is spread 
evenly on surfaces of constant Hq. 

We start by decomposing the true distribution function of par¬ 
ticles in phase space, /, in terms of a spherical background /o and 
a perturbation Sf. To make sure the split between spherical back¬ 
ground and aspherical perturbation is uniquely defined, we take fo 
as the distribution function obtained when we perform a naive anal¬ 
ysis averaging out the aspherical contribution: 

M J) =j^fd 3 @f(J,@). (9) 

By Jean’s Theorem, fo is an equilibrium distribution function in the 
spherical background because it is constructed from spherical in¬ 
variants J alone. Analogous to equations 0 and 0 one can write 
the full distribution function / as 


/(J,0)=/o(J) + 5/(J,0) 

= MJ) + '£8f n (J)e in& . (10) 

n 


The whole / is to be in equilibrium in the true system, 
df/dt = 0. We can turn this into an explicit condition on / using 
the the collisionless Boltzmann equation. 


0 = 


df 

dt 


\h f i = — d L_ dl L ^j_ 
[ ~ d® dj dj'd® 


(ii) 


Expanding to linear order in H and / gives the condition 


£ (flo (J)-nSf n (J) - ^j-nSH n (J^j e in& = 0. (12) 


The different 0 dependence of each term in the sum means that 
the term in brackets must be zero for each n. In particular, for the 
resonant terms where $~2y ■ = 0 one has the condition 

SH rl± (J) = 0. (13) 


A sufficient condition for stability of Sf is therefore that fo is a 
function only of Ho, since then dfo/dJ = J"2o d/o/d//o and conse¬ 
quently the dot product in equation {T3} vanishes. 

This is the core result claimed at the start of the section: fo = 
fo(Ho), i.e. the distribution function implied by a spherical analysis 
appears to be ergodic. It is not a necessary condition for achieving 
equilibrium, since for any given aspherical system certain SH n± 
will be zero. Rather, the result should be read as applying to the 
maximally stable distribution function - a distribution function that 
does not evolve under any linear perturbation to its potential. 
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4 A. Pontzen et al 


We again emphasise that the distribution function /o is fic¬ 
tional. There is no sense in which the true distribution function, /, is 
actually ergodic. The statement is about how the system appears to 
be when it is analysed using spherically averaged quantities, equa¬ 
tion & Yet, it establishes a way in which we can understand these 
spherically-averaged quantities in a systematic way - the system is 
most stable if it appears ergodic, regardless of what the underlying 
dynamics is really up to. In the remainder of this paper we will refer 
to such systems as ‘spherically ergodic’. 


2.3 Testable predictions 


We have established that systems which appear to be ergodic in a 
spherical analysis are maximally stable. Now we need to devise a 
connection to observable or numerically-measurable quantities. 

A distribution function /(Ho) that is truly a function of energy 
alone has an isotropic velocity distribution ( |Binney & Tremaine| 
2008 1 . To test for isotropy, one calculates /3 (r) according to the 
usual spherically-averaged definition 


P(r)=\ 


(*?)(>•) 

2 (v?)(r) 


(14) 


where v r is a particle’s radial velocity, v t its tangential velocity and 
the angle-bracket averages are taken in radial shells. For a pop¬ 
ulation on radial orbits, /3(r) = 1; conversely for purely circular 
motion, j8 (r) = — °°. Between these two extremes, an ergodic pop¬ 
ulation has j3 (r) = 0 i Binney & Tremaine|2008f . 

Intuitively, a spherically ergodic system (in the sense de¬ 
fined in the previous section) should therefore be approximately 
isotropic. However one has to handle that expectation with a little 
care because the true population / is not ergodic and the measured 
velocity dispersions, even in spherical polar coordinates, may in¬ 
herit information from / that is not present in /q. 

Instead we will now construct a more rigorously justifiable, 
slightly different statement that still connects spherically ergodic 
populations to velocity isotropy. Measuring the mean of any func¬ 
tion of the spherical actions q(J), we obtain 

f d 3 J d 3 0/( J ,&)q(J) = j d 3 _// 0 (H 0 (J)) 9 (J), (15) 

an exact result. Therefore any statement about averages over spher¬ 
ical actions automatically knows only about /o - the spherical part 
of the distribution function. This allows us to derive unambiguous 
implications of a spherically ergodic population. 

The most familiar action is the specific scalar angular momen¬ 
tum j. Because it is a scalar for each particle, averages over this 
quantity do not express anything about a net spin of the halo but 
rather about the mix of circular and radial orbits, just like the tradi¬ 
tional velocity anisotropy. Radially-biased populations have (j) ~ 0 
whereas populations on circular orbits have (j) = j c , where j c is 
the maximum angular momentum available at a given energy. So, 
velocity anisotropy can be conveniently represented in terms of the 
mean scalar angular momentum. 

We can go further and calculate a function, (j) ( E ), where the 
average is taken only over particles at a particular specific energy. 
This quantity can be represented in terms of the ratio of two inte¬ 
grals of the form m 

JJJdJ r djdj z MH 0 )jS(H 0 -E) 


(j)(E) = 


ffJdJ r djdj z MH 0 )8(H 0 -E) ' 


(16) 


The triple integral ranges over the physical phase space coordi¬ 
nates: 0 Z /- < °°, 0 Z 7 < °°, —j Z j- Z j. One can immediately 


perform the j- integrals; then the J r integral can be completed by 
changing variables to Hq (recalling <9Ho/ dj r = £ 2 ,) and consuming 
the 8 function. After this manipulation j can only range between 
0 and j c (E) where j c (E) is the specific angular momentum cor¬ 
responding to a circular orbit with specific energy E ; there are no 
physical orbits with more angular momentum at the specified E. 
The final, exact result is: 

fjc( E ) . 0 / rj c (E) 

(j)(E)= l o d j£l r (E,j)~ l j 2 /J o djn r (E,j)~ l j. (17) 

We now have a firm prediction for spherically ergodic populations. 
Namely, if we bin particles in E and measure ( j)(E ) in each bin, 
the results should be predicted by equation ITT] which is a function 
only of the potential (through £l r ). Equation (| 1 7|> does not exhaust 
the possible tests for spherical ergodicity, but it is sufficient for our 
present exploratory purposes. 

For smooth potentials, £ 2, •(£,_/') varies very little between 
j = 0 and j = j c , and one can approximate it very well as a func¬ 
tion of E alone (Appendix [A3}. In this case, the integrals follow 
analytically and one has the result 

(j)(E)~lj c (E). (18) 

This is a helpful simplification to set expectations, but throughout 
this paper when showing the spherically ergodic limit, we will use 
the exact expression given by equation G3- 


3 EXAMPLE APPLICATIONS 

So far we have motivated and derived a formal result that aspheri- 
cal systems are most stable when they appear ergodic in spherical 
coordinates. We derived one practical consequence for the angular 
momentum distribution, equation (17) . which in an approximate 
sense states that the velocity distribution will appear isotropic. We 
are now in a position to test whether numerical simulations actually 
tend towards this maximally-stable limit in a variety of situations, 
beginning with cosmological collisionless dark matter halos. 


3.1 Cosmological dark matter halos 

Let us re-examine the three high-resolution, dark-matter-only zoom 
cosmological simulations used in the analysis of |Pontzen & Gow| 
|emato| (2013) . The three each have several million particles in 
their z = 0 halos which correspond in turn to a dwarf irregular, 
L* galaxy and cluster. The force softening e, virial radius ri oo (at 
which the mean density enclosed is 200 times the critical density) 
and virial masses M 200 are 65. 170, 690pc; 98, 301, 1430kpc and 
2.8 x 10 10 , 8.0 x 10 11 , 8.7 x 10 13 M Q respectively. For further de¬ 
tails see |Pontzen & Governato| l f2013| >. 

The upper panel of Figure|2]shows the anisotropy for our cos¬ 
mological halos. To compare the three directly, we scale the ra¬ 
dius by r max (respectively 27, 57 and 340kpc), the radius at which 
the circular velocity (GM(< r)/r) l P reaches its maximum, v ma x 
(= 56, 150 and 610kms _1 ). We restrict attention to the region well 
within the virial radius; here, the anisotropy /3 (r) typically lies be¬ 
tween the purely radial and isotropic cases (e.g. Bellovary et al.| 
[20081 |Navarro et al,|2010) . 

We now want to link this relatively familiar velocity 
anisotropy to the alternative ( j)(E ) statistic that was directly pre¬ 
dicted by the spherically ergodic property in Section [2] For each 
particle we calculate the specific energy E = ar/2 + <I> (|a;|), where 
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r/fmax 

0.0 0.5 1.0 1.5 2.0 



E/Vmax 

Figure 2. The velocity anisotropy of the inner parts of three sample high- 
resolution cosmological dark matter halos (simulated without baryons), 
plotted as a function of radius (upper panel) and in energy shells (lower 
panel). The upper panel shows the classic velocity anisotropy /3(r) defined 
in the text, for which a purely radial population has /3 = 1 and a population 
on circular orbits /3 = — °o. The lower panel shows our alternative in energy 
space which can be more precisely related to the theoretical arguments pre¬ 
sented in Section[2j here ( j){E)/j c \ rc = 0 for radial orbits and 1 for circular 
orbits, and approximately 2/3 for an isotropic population. Both panels show 
that the halos have near-isotropic orbits with a slight radial bias. The range 
of the two plots is roughly comparable, but we caution that the mapping 
from r to E is not unique (see Figure [3}. 


x is the vector displacement from the halo centre. To make this 
quantity agree exactly with Ho in the terminology of Section [2] we 
ignore asphericity when calculating the potential energy, defining 
it as 

„ 9) 

Jo r z 

where M(< /) is the mass enclosed inside a sphere of radius /. 
(The numerical integration is performed by binning particles in 
shells of fixed width e, chosen to coincide with the force soften¬ 
ing in the simulation; within these bins the density is taken to be 
constant.) The physics is invariant if a constant is added to the po¬ 
tential; we have chosen to fix its scale by setting <5(0) = 0. 

For each particle we also calculate the specific angular mo¬ 
mentum j = \x x x\, and the specific angular momentum of a cir¬ 
cular orbit at the same energy, j c kc(E) which is given by simulta¬ 
neously solving 

•2 '2 

E = <f>(r) + J -Sf and (20) 

2 r l 

to eliminate r in favour of y c ; rc . 

We plot (j)(E)/jarc(E) in bins containing 1000 particles 
each in the lower panel of Figure [2] To facilitate comparison with 
the top panel, a population on purely radial orbits would have j8 = 1 
and (j) /y'circ = 0, whereas a purely circular distribution function 
corresponds to /3 = — °o or (j) /y c ; rc = 1. Isotropic, purely spherical 
populations have /5(r) =0 and (y)/y'circ — 2/3 as discussed at the 
end of Section |23l When compared against each other in this way, 
the two panels agree well: the populations are on near-isotropic or¬ 
bits with a slight radial bias. 

Quantitatively, how well do these results agree? As a guide¬ 
line, we can compare results for models with constant anisotropy. 



Figure 3. The relationship between r/r max and T/V/ nax for circular orbits 
(dashed line), radial orbits at apocentre (solid line) and particles in the 
‘MW’ run (density shows the number of particles with energy E at each ra¬ 
dius r). The mapping between energy and radius is fuzzy, so that anisotropy 
at high E can easily contaminate j8 (r) at small r. 


From |Binney & Tremaine |2008| , a distribution function f(j,E) = 
j-^m generates a constant anisotropy /5(r) = j8. We can cal¬ 
culate the connection to the new statistic by generalising the rea¬ 
soning of Section[T3| with the result that 


(ME) 


jj^djQ.- l j 2 - 2 P 

Jo drc dyTV 1 y 1 ~ 2 /* 


2 — 2/5 . 


( 21 ) 


where the first result is exact and the second follows from assum¬ 
ing £l r is independent of j (which is an excellent approximation; 
see Appendix |A3| >. Consistent with equation ( |18| >, [5=0 gives 
(y) (E)/ y c j rc ~ 2/3. But as the system becomes more radially biased 
we can now calculate that, for example, [5 = 0.2 corresponds to 
(j) (E) /y c irc — 0.62. Despite the various approximations involved, 
these values therefore correctly relate the values of /5 in the top 
panel of Figure[2]with the (y) results in the lower panel. 

That said, a detailed comparison as a function of radius is hard 
because particles at a given radius r have a wide spread of ener¬ 
gies E. Figure[3]illustrates this relationship for the ‘MW’ halo. The 
density shows the probability of a particle at radius r also having 
specific energy E , p(E\r). The minimum E at each r is set by <I>(r), 
which gives the energy of a particle at apocentre (solid line). A 
more typical E is given by the energy of a circular orbit at r (dashed 
line), and this gives some intuition for mapping results from the top 
panel of Figure [2] onto the bottom panel. However, any E exceed¬ 
ing 4>(r) is theoretically possible. So j6 at any given radius actually 
represents an average over particles of many different energies. 


3.2 The classical radial orbit instability 

Having established that, loosely speaking, ( j)(E ) represents the 
anisotropy in energy shells in the same way that [5 ( r ) does in radial 
shells, we can return to our prediction (Section [23} for the former 
quantity, which is shown by the dashed line in the lower panel of 
Figure [2] The prediction is almost, but not quite, satisfied in cos- 
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r/fmax 


0.0 r, 0.5 1.0 1.5 2.0 



Figure 4. The p(r) velocity anisotropy (upper panel) and its equivalent 
in energy space ( j)(E ) (lower panel, as Figure [2j for the evolution over 
time of a halo that is initially in equilibrium, but unstable against the radial 
orbit instability. The initial conditions are represented by the dark blue line, 
t = 0; there is a delay of 12/ ( j VII before any significant evolution can be 
seen (light blue line). Subsequent lines are shown every 4/<| yn . Once the 
instability kicks in and generates aspherical perturbations, there is a rapid 
evolution towards the predicted spherically ergodic limit (lower panel). 


mological dark matter halos. Since the condition is only reached in 
a maximally-stable object, approximate agreement is an acceptable 
situation. 

Because cosmological halos initially form from near-cold col¬ 
lapse, the radial orbit instability (van Albada 1982; Barnes et al. 
|1986| |Saha|199l| |Weinberg|| 1991 | is invoked to explain how the 

radially infalling material gets scattered onto a wider variety of or- 
bits|MacMillan et HTTj ( |2()()6^ ; | Bellovary et al.|f2008| l; [Barnes et al.] 
(2009), isotropising the velocity dispersion. Most tellingly, numer¬ 
ical experiments by |Huss et ah] < | 1999|) ; |MacMillan et al.| j2006| 
show that the suppressing the instability (by switching off non- 
radial forces) results in a qualitatively different density profile as 
the end-point of collapse. We will now show that the isotropisation 
of velocity dispersion during the radial orbit instability can be in¬ 
terpreted in terms of an evolution towards stability in the terms of 
Section[2] 

Consider what happens to a halo that is intentionally designed 
to be unstable. First, we will show the classic radial orbit instability 
(ROI) at work by constructing a spherical halo with particles that 
are on radially-biased orbits. We initialise our particles such that 
they solve the Boltzmann equation and so are stable in exact spher¬ 
ical symmetry. In practice, however, the strong radial bias means 
that any slight numerical noise will trigger the ROI. By initialis¬ 
ing an unstable equilibrium in this way, we avoid confusion from 
violent relaxation processes associated with out-of-equilibrium col¬ 
lapse ( |Lynden-Bell| 1 967| l . 

Specifically, the initial conditions are set up in a similar fash¬ 
ion to|Read et al.[ (2006 ), with particle positions drawn from a gen¬ 
eralised Hemquist density profile (Hernquist 1990, Dehnen 1993): 


p{r) = 


Af(3 — y) fr\~Y 
4 7ta 3 


1 + ■ 


7—4 


( 22 ) 


which has a circular velocity reaching a peak at r max = (2 — y)a 
and implies the gravitational potential 


®w=^,[ (i +«Ar 2 - 1 ]. (23) 


We choose y = 1 to roughly mimic an NFW halo (Navarro et al. 
|l997j > in the innermost parts of interest. The velocities are sampled 
(using an accept-reject algorithm) from a numerically-calculated 
Osipkov-Merritt distribution function (Binney & Tremaine |2008| l: 

m [(i+'W 2 m)pm«>)] (24) 

where the parameter Q is defined by 


Q = E+ k- 


(25) 


The value of r a is known as the anisotropy radius because the ve¬ 
locity anisotropy is given by 


P( r ) — 


1 


l + rl/r 2 ’ 


(26) 


showing that jS (r) ~ 0 (isotropic) for r <C r a and ]3 (r) ~ 1 (radial) 
for r>r a . We used the minimum value of r a for which the distribu¬ 
tion function is everywhere positive, making the orbits as radially 
biased as possible; for 7=1, this is r a ~ 0.21 ( |Meza & Zamorano] 
1997). We draw 10® particles and evolve the system using RAM¬ 
SES (Teyssier |2002} , with mesh refinement based on the number of 
particles per cell, resulting in a naturally adaptive force softening 
reaching a minimum of e = 90pc. 

Our expectation that numerical noise triggers the ROI is borne 
out by the numerical experiment. The upper panel of Figure [4] 
shows the radial anisotropy /3 (r) over time. We have defined a 
single dynamical time, t ( j yn , at the peak of the velocity curve so 
that t ( j yn = r max /v max . The six solid lines show the population at 
t = 0, 12,16,20,24 and 28 dynamical times. At first, /3 (r) appears 
stable, but suddenly after 16 ?d yn it becomes considerably more 
isotropic. Over the same timescale, asphericity in the potential de¬ 
velops. To demonstrate this, we determine the inertia tensor of the 
entire density distribution and calculate the ratio of the principal 
axes; at t = 0, the ratio is 1.0 by construction. By 16tdyn th e ratio 
is 1.8. It stabilises at around 26?d yn with a ratio of 4.3. 

This is symptomatic of the classic radial orbit instability in 
action. We can follow the same process from our energy/angular- 
momentum standpoint in the bottom panel of Figure [4] The initial 
conditions (At = 0) have a kink in them, with most regions of en¬ 
ergy space appearing radially biased but some showing a slight cir¬ 
cular preference (around E ~ 1.5ty nax ). We verified that this is an 
artefact of the Opsikov-Merritt construction and is consistent with 
the radial /3(r) being radially biased everywhere; recall Figure [3] 
shows that the relationship between r and E is non-trivial. 

As time progresses, the (j) (E) distribution correctly tends to¬ 
wards the spherically ergodic (SE) limit, as predicted. The SE limit 
is attained to good accuracy for energies E < by At = 20rd yn . 
At larger energies, it is likely difficult to achieve SE because of the 
long time-scales and weak gravitational fields involved. Like (j), 
jS (r) at 28?dyn shows an isotropic distribution at the centre (r —¥ 0); 
but /3 (r) > 0 everywhere for r > 0.1 r max . We verified that this con¬ 
tinuing radial bias throughout the halo is produced by a number 
of high-energy, loosely-bound particles plunging through, i.e. the 
high-E particles from the lower panel in Figure |4]pervade all radii 
in the upper panel. 

Let us briefly recap: what we have so far is a new view 
of an existing phenomenon. We have recast the ROI and its im¬ 
pact on spherically-averaged quantities as an evolution towards an 
analytically-derived maximally stable class of distribution func¬ 
tions. Viewed in energy shells instead of radial shells, the distinc¬ 
tion between the regions that reach the SE limit and the slowly- 
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Figure 5. The mean angular momentum of particles as a function of energy, 
like the lower panel of Figure [4] — but now for a subpopulation in a dark 
matter halo that is globally in equilibrium. At At = 0, we select particles on 
predominantly radial (j/ / CMC < 0.2) orbits; at later times, the subpopulation 
mean evolves back towards the population mean. The lines from bottom 
to top show the state at selection and after 0.2, 0.4, 0.6, 1.0 and 4.0f,j yn 
respectively. 

evolving, loosely-bound regions that remain radially biased is con¬ 
siderably clearer. Now, because our underlying analytic result is 
derived without requiring the population to be self-gravitating, we 
can now go further and consider a different version of the radial 
orbit instability - and show that it continues to operate even when 
a system is in a completely stable equilibrium. 

3.3 The continuous radial orbit instability 

Our next target for investigation is to broaden the conditions: ac¬ 
cording to the results of Section [2] a subpopulation of particles 
should undergo something much like the radial orbit instability 
even when the global potential is completely stable. A suitable 
name for this phenomenon would seem to be the “continuous ra¬ 
dial orbit instability”, since it continues indefinitely after the global 
potential has stabilised. 

We perform another numerical experiment to demonstrate the 
effect. First, to avoid confusion from cosmological infall and tidal 
fields, we create a stable, isolated, triaxial cosmological halo by 
extracting from our cosmological run a region of 3r2oo around our 
‘Dwarf’. We then evolve this isolated region for 2 Gyr to allow any 
edge effects to die away, and verify that the density profile out to 
7*200 is completely stable. As before we define a dynamical time for 
the system of r dyn = r max /v max = 470Myr. 

After the 2 Gyr ~ 4t dyn has elapsed, we select all particles with 
j < 0.2 j QiIc . These particles are, at the moment of selection, on 
preferentially radial orbits. We trace our particles forward through 
time, measuring ( j)(E ) in each snapshot. The results are shown in 
Figure [5j for various times between At = 0.0 and 4.0/ dyn after se¬ 
lection. Over this period, the mean angular momentum significantly 
increases towards the spherically ergodic limit at every energy. The 
changes are much faster at low energies where the particles are 
tightly bound and the local dynamical time is short compared to 
the globally-defined f dyn . 

The evolution is rapid until 4t dyn after which the ( j)(E ) re¬ 
mains near-static (except at E > 2.5v^ ax where it continues to 


r/r max 

0.0 0.1 0.2 0.3 0.4 0.5 



E/Vmax 

Figure 6. The j8(r) (upper panel) and ( j)(E ) (lower panel) relations for the 
stellar populations of dwarf satellites around a Milky-Way-like central. In 
the upper panel arrows indicate the location of , the radius enclosing half 
the stellar mass. The profile is plotted exterior to 3e where e is the softening 
length. Because of the tidal interaction with the parent galaxy, the outer 
parts of each object are out of equilibrium and display biases ranging from 
strongly radial to circular. The lower panel shows the same populations in 
energy space. At low energies, tightly bound particles are now seen to be 
close to the spherically ergodic limit. At higher energies, a spread is seen 
but not as large as would naively be expected from the (r) relation. 

slowly rise). Eventually the subpopulation has ( j)(E ) ~ 0.5y’circCE) 
independent of E. This establishes that particles at low angular mo¬ 
mentum within a completely stable, unevolving halo, automatically 
evolve towards higher angular momentum. After a few dynamical 
times, their angular momentum becomes comparable to that of a 
randomly-selected particle from the full population, although with 
a continuing slight radial bias. 

One can understand this incompleteness in a number of equiv¬ 
alent ways. Within our formal picture, it arises from the fact that 
only certain 8H n are non-zero in equation ( [ 1 2) . More intuitively, 
the continuing process of subpopulation evolution is being driven 
by particles on box orbits that change their angular momentum at 
near-fixed energy - but not all particles are on such orbits, and so 
some memory of the initial selection persists. 

3.4 Observations of dwarf spheroidal galaxies 

In the previous section, we established that subpopulations on ini¬ 
tially radially-biased orbits evolve towards velocity isotropy even 
if the global potential is stable. We explained this in terms of our 
earlier calculations, and now turn to how the results might impact 
on observations. 

An important current astrophysical question is how dark mat¬ 
ter is distributed within low-mass galaxies: this can discriminate 
between different particle physics scenarios jPontzen & Governato| 
|2014| >. The smallest known galaxies, dwarf spheroidals surround¬ 
ing our Milky Way, in principle provide a unique laboratory from 
this perspective. But determining the distribution of dark matter is 
a degenerate problem because of the unknown transverse velocities 
of the stellar component (e.g. [Charbonnier et aLj|2011| l. Further¬ 
more the use of spherical analyses seems inappropriate since the 
underlying systems are known to be triaxial (e.g. |Bonnivard et al.| 

[20lSl >. 
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This is exactly the kind of situation we set up in our initial cal¬ 
culations (Section[2]l: an aspherical system being analysed in spher¬ 
ical symmetry. We have shown that the results apply both to tracer 
and self-gravitating populations. So, we can go ahead and apply it 
to the stars in observed systems: the stellar population of a dwarf 
spheroidal system will be maximally stable if it appears ergodic 
in a spherical analysis. This implies a strong prior on what spher¬ 
ical distribution functions are actually acceptable and therefore, in 
principle, lessens the anisotropy degeneracy. 

This idea warrants exploration in a separate paper; here we 
will briefly test whether the idea is feasible by inspecting some 
simulated dwarf spheroidal satellite systems. In particular we use 
a gas-dynamical simulation of “MW” region (see above) using ex¬ 
actly the same resolution and physics as |Zolotov et al. | ( [20 1 2ft ; see 
that work for technical details (although note the actual simulation 
box is a different realisation). 

We analyse all satellites with more than 10 4 star particles, 
which gives us objects lying in the ranges 5.4 x 10 7 < M* < 
4.9 x 1O 8 M 0 , 28 < v max < 60kms -1 and 4.6 < r max < 9.6kpc. 
We first verified that these do not host rotationally supported disks; 
we then calculated p(r) profiles between 3e ~ 0.5kpc and 0.6r max . 
The stellar half-light radius is 0.15 < D/ 2 /V max < 0.3, so our cal¬ 
culations extend into the outer edges of each visible system. 

The results are shown in the upper panel of Figure [6] with ar¬ 
rows indicating rq The outer regions display a range of different 
behaviours from strongly radial to circular. However, in each case 
the P ( r ) is nearly isotropic as r —> 0. This is consistent with a view 
in which the centres have achieved stability while the outer parts 
are being harassed by the parent tides and stripping. 

The picture is reinforced by considering the { j)(E ) statistic 
(lower panel of Figure[6]l. Here, stars on tightly bound orbits (small 
E) are very close to the spherically ergodic limit. In the less-bound 
regions, the agreement worsens. However a quantitative analysis 
using the approximate relation © shows that the ( j)(E) relation 
stays much closer to the spherically ergodic limit than the p(r) 
relation naively implies. This is probably because the shape of p (r) 
is in part determined by out-of-equilibrium processes, coupled to 
the multivalued relationship between r and E (Figure [3}. 

In any case these results suggest that we should adopt priors 
on spherical Jeans or Schwarzschild analyses that strongly favour 
near-isotropy in the centre of spheroidal systems - or, more accu¬ 
rately, strongly favour near-ergodicity for tightly bound stars. In 
future work we will expand on these ideas and apply them to ob¬ 
servational data, since they could lead to a substantial weakening 
of the problematic density-anisotropy degeneracies. 

3.5 Dark matter cusp - core transitions 

In jPontzen & Governato| (20 12[ henceforth PG12) we established 
that dark matter can be redistributed when intense, short, repeated 
bursts of star formation repeatedly clear the central regions of a 
forming dwarf galaxy of dense gas (see also jRead & Gilmore |2005[ 
[Mashchenko et al.|2006| l. The resulting time-changing potential im¬ 
parts net energy to the dark matter in accordance with an impulsive 
analytic approximation. This type of activity has now been seen 
or mimicked in a large number of simulations, allowing glimpses 
of the dependency of the process on galactic mass, feedback type 


and efficiency 

Governato et al. 2012; Penarrubia et al. 2012; Mar- 

tizzi et al. 2012 

Teyssier et al. 2013 

ZolotovetalJ201^; Garrison- 

Kimmel et al. 2013; Di Cintio et al. 

2014; Onorbe et al. 2015 ). 


We based our PG12 analysis on 3D zoom cosmological simu¬ 
lations, but our analytic model assumed exact spherical symmetry. 


r/kpc 




E/km 2 s " 2 

Figure 7. The velocity anisotropy (/3 (r), upper panel) and angular momen¬ 
tum (( j)(E ), lower panel) for cusped and cored halos from PG12; note the 
much-expandedy-axis scales relative to previous figures, which are required 
to highlight the differences. The cored cases (solid lines) are almost per¬ 
fectly spherically ergodic (lower panel) and hence isotropic (upper panel). 
This contrasts with the cusped case (dotted lines) which has a slight but sig¬ 
nificant radial bias (seen as high P in the upper panel and low j in the lower 
panel). 

By definition, in the analytic model all particles conserve their an¬ 
gular momentum at all times. Taken at face value, energy gains 
coupled to constant angular momentum would leave a radially bi¬ 
ased population in the centre of the halo. 

The top panel of Figure [7] shows the measured velocity 
anisotropy in the inner 5 kpc for the zoom simulations in PG12 at 
z = 0; the solid line shows the feedback simulation which has de¬ 
veloped a core, whereas the dotted line shows the dark-matter-only 
simulation which maintains its cusp. The difference between the 
two cases is the opposite to that naively expected: the centre of the 
cored halo has a more isotropic velocity dispersion than the centre 
of the cusped halo. The lower panel shows the equivalent picture in 
energy space (noting that <f>(0) = 0 and <F(5kpc) = 3790km 2 s~ 2 ). 
For clarity only a small fraction of the 0 < j/ j c ; rc < 1 interval is 
now plotted; the differences between the cusped (dotted) and cored 
(solid) lines are relatively small, but significant. At all energies the 
cored profile has more angular momentum than the cusped profile; 
in fact, it lies right on the spherically-ergodic limit (dashed line) 
whereas the cusped profile is biased by around 0.06 to radial or- 
biti[]] 

While the PG12 model correctly describes the energy shifts, 
it misses the angular momentum aspect of core creation, which 
has been emphasised elsewhere (e.g. Tonini et al .j2006| ). We will 
now show that a complete description of cusp-core transitions in¬ 
volves two components: energy gains consistent with PG12, and 
re-stabilisation of the distribution function consistent with the de¬ 
scription of Sections [2.2| and [3~3] We again use the RAMSES code 
to follow isolated halos with particle mass 2 x 1O 4 M 0 ; however we 
adapted the code to add an external, time-varying potential to the 
self-gravity. 

1 The dashed line for the isotropic j/j civc in the lower panel is calculated 
using the cored potential. However the dependence on potential is very 
weak, as discussed in Section [5] calculating it instead with the cusped po¬ 
tential leads to differences of less than one percent at every energy. 
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Figure 8. The PG12 mechanism is reproduced by adding an external poten¬ 
tial (to mimic ‘gas’) to a DM-only simulation for the time intervals indicated 
by the grey bands. We measure the response of a completely spherical halo 
(dotted lines) and a realistically triaxial halo (solid lines) as described in the 
text. In both cases we track the mean energy (top panel) and angular mo¬ 
mentum (middle panel) of the 0.1% of particles that start out most bound. 
We also measure the log slope of the density profile at each time output 
(bottom panel). The energy shift (top) predicted by PG12 is found to hold 
in both cases. In the spherical case, the angular momentum remains constant 
up to discreteness effects; it therefore drops relative to the circular angular 
momentum (middle panel). Conversely, in the triaxial case, the continuous 
radial orbit instability causes the angular momentum to rise in proportion 
to jcirc- Only when the angular momentum is allowed to rise does a sig¬ 
nificant core form, meaning that the spherical case unphysically suppresses 
core formation (lower panel). 

We took our stabilised, isolated, triaxial dark matter halo ex¬ 
tracted as described in Section [33| and created a sphericalised ver¬ 
sion of it as follows. For each particle we generate a random ro¬ 
tation matrix following an algorithm given by |Kirk| (JT992J, then 
multiply the velocity and position vector by this matrix. Finally, we 
verified that the final particles are distributed evenly in solid an¬ 
gle, and that the spherically-averaged density profile and velocity 
anisotropy is unchanged. The triaxial and spherical halos are both 
NFW-like and stable over more than 2 Gyr when no external poten¬ 
tial is applied. 

For our science runs, we impose an external potential corre¬ 
sponding to 10 8 Mq gas in a spherical ball lkpc in radius, dis¬ 
tributed following p r~ 2 ; this implies a potential perturbation 
of A<t> = —700km 2 s~ 2 at 500pc, for instance. The potential in¬ 
stantaneously switches off at lOOMyr, back on at 200Myr, off at 
300Myr and so forth until it has accomplished four “bursts”. The 
period and the mass in gas is motivated by Figures 1 and 2 of PG12 
and Figure 7 of T+13. We also tried imposing potentials with dif¬ 
ferent regular periods and with random fluctuations, none of which 
altered the behaviour described below. 

The top two panels of Figure [8] show the time-dependent be¬ 
haviour of the central, most-bound 0.1 % of particles. The triaxial 
and spherical halo results are illustrated by solid and dashed lines 
respectively. The top panel shows how the coupling of the external 
potential is similar; in particular, it results in a mean increase in spe¬ 
cific energy of A E ~ 200 km 2 s -2 for particles in both cases. The 
final shift in the spherical case is very slightly larger than that in the 
triaxial case, a difference which is unimportant for what follows (it 
would tend to create a larger core if anything). 


The middle panel displays the mean specific angular momen¬ 
tum j of the same particles as a fraction of the circular angular 
momentum j c i TC (E). In the spherical case (dashed line), this quan¬ 
tity drops significantly because j is fixed but j c „. c is rising over time 
(since it increases with E). By contrast in the triaxial case j / j c [ rc 
returns to its original value, meaning that j has risen by the same 
fraction as j C uc(E). This is dynamic confirmation of the discussion 
earlier in this Section: the spherically symmetric approach predicts 
an increasing bias to radial orbits - whereas in the realistic aspher- 
ical case, the stability requirements quickly erase this bias. 

The bottom panel of Figure [8] shows the measured density 
slope a = dlnp/dlnr at 500pc for the two experiments. Both 
start at a ~ —1.5; the triaxial case correctly develops a core (with 
a ~ —0.1) whereas the completely spherical case maintains a cusp 
(with a ~ —0.9). This discrepancy is causally connected to the an¬ 
gular momentum behaviour: the mean radius of a particle increases 
with increasing angular momentum j, even at fixed energy E, so a 
typical particle migrates further outwards in the triaxial case com¬ 
pared to the spherical. 

We can therefore conclude that asphericity is a pre-requisite 
for efficient cusp-core transitions. Flowever there is a subsidiary 
issue worth mentioning. After about 3 Gyr we find that the spheri¬ 
cal halo autonomously does start increasing (j) / j c ; rc and the dark 
matter density slope becomes shallower. This is because the poten¬ 
tial fluctuations have generated a radially biased population which 
is unstable, and a global radial orbit instability is triggered by nu¬ 
merical noise over a sufficiently long period (as in Section [3~2] i- 

In fact, with ‘live’ baryons rather than imposed fluctua¬ 
tions, there are aspherical perturbations which accelerate the re¬ 
equilibration process further and renew a global radial orbit in¬ 
stability (Section [3 .2[ >, encouraging the entire population towards 
spherical ergodicity. As we saw in Figure [7] the cored halo from 
PG12 has an almost perfectly spherically-ergodic population for 
E < 3 000km 2 s -2 (unlike the cusped case). We verified that this is 
also true of T+13. Generating a core through potential fluctuations 
seems to complete a relaxation process that otherwise freezes out 
at an incomplete stage during collisionless collapse. 


4 CONCLUSIONS 

Let us return to the original question: how much do inaccuracies 
inherent in the spherical approximation really matter in practical 
situations? The answer is, unfortunately, that “it depends”; but we 
can now distinguish two cases as a rule of thumb: 

1) For dynamical calculations or simulations, the inaccuracy mat¬ 
ters a great deal. Neglecting the aspherical part of the potential 
unphysically freezes out the radial orbit instability and related 
effects, so can lead to qualitatively incorrect behaviour. 

2) For the analysis of observations or simulations in equilibrium, 
the assumption is far more benign - it is, in fact, extremely pow¬ 
erful when handled with care. The underlying aspherical system 
and the fictional spherical system both appear to be in equilib¬ 
rium; the mapping between the two views yields striking insight 
into (for example) spheroidal stellar distribution functions and 
dark matter halo equilibria. 

These conclusions are based on the fact that, when aspherical 
systems are analysed in spherical coordinates, there is an attrac¬ 
tor solution for the spherically-averaged distribution function /o - 
namely, it tends towards being ergodic (i.e. /o is well-approximated 
as a function of energy alone). We demonstrated this using Flamil- 
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tonian perturbation theory (see Section [2] and Appendix |Bj, and 
subsequently used the term “spherically ergodic” (SE) to describe 
a distribution function / with the property that its spherical average 
/o is ergodic in this way. 

The result follows because the orbits do not respect spheri¬ 
cal integrals-of-motion such as angular momentum. Note, however, 
that the physical orbits may still possess invariants in a more ap¬ 
propriate set of coordinates. The apparent chaotic behaviour and 
tendency of particles to spread evenly at each energy is a helpful 
illusion caused by adopting coordinates that are only partially ap¬ 
propriate. 

In Section [3] we showed that the idea has significant explana¬ 
tory power. First, we inspected the selection of equilibrium in tri- 
axial dark matter halos (Section [3. 1 [ 1 which led us to consider the 
classical radial orbit instability (Section[32]). We demonstrated that 
the instability naturally terminates very near our SE limit (lower 
panel, Figure [ 4 ). Particles at high energies have long dynamical 
times which causes them to freeze out: they evolve towards, but do 
not reach, the SE limit. Because these high-energy particles often 
stray into the innermost regions (Figure|3]l, the velocity anisotropy 
/3 (r) continues to display a significant radial bias at all r after the 
instability has frozen out. Moreover in the case of self-consistently 
formed cosmological halos (Figure|2]l, even at low energies there is 
a slight radial bias. The bias is only erased when a suitable exter¬ 
nal potential is applied ( e.g. during baryonic cusp-core transforma¬ 
tions), forcing the system to stabilise itself against a wider class of 
perturbations than it can self-consistently generate (Figure |7J; we 
will return to this issue momentarily. 

One novel aspect of our analysis compared to previous treat¬ 
ments of the radial orbit instability is that it applies as much to 
tracer particles as to self-gravitating populations. As a first exam¬ 
ple, in Section [33] we demonstrated how a subset of dark matter 
particles chosen to be on radially-biased orbits mix back into the 
population (Figure [5}. This is the case even for halos that, as a 
whole, are in stable equilibrium; therefore we referred to the phe¬ 
nomenon as a ‘continuous’ radial orbit instability. 

The same argument implies that stars undergo the radial orbit 
instability as easily as dark matter particles. In particular, without 
a stable disk to enforce extra invariants, dwarf spheroidal galaxies 
likely have stellar populations with a near-SE distribution (Section 
|3.4| Figure [6]». This provides a footing on which to base spherical 
Jeans or Schwarzschild analyses of observed systems: it implies an 
extra prior which can be formulated loosely as stating that j8 (r) —>■ 0 
as r —¥ 0. Such a prior could be powerful in breaking the degen¬ 
eracy between density estimates and anisotropy (e.g. |Charbonnier| 
|et al .|20 In turn tightening limits on the particle physics of dark 
matter [Pontzen & Governato|2014) . 

As a final application, we turned to the question of the bary¬ 
onic processes that convert a dark matter cusp into a core (Section 
|3.5[ >. Angular momentum is gained by individual particles during 
the cusp-flattening process ( fTonini et al.|2006; > but our earlier work 
(especially PG12) has focussed primarily on the energy gains in¬ 
stead. In Figure [8] we see that, in spherical or in triaxial dark mat¬ 
ter halos, time-dependent perturbations (corresponding to the be¬ 
haviour of gas in the presence of bursty star formation feedback) 
always lead to a rise in energy. However, in the spherical case the 
mean angular momentum as a fraction of j ckc drops because the an¬ 
gular momentum of each particle is fixed, whereas j ckc rises with 
E. Only in the triaxial case does (j) /y c ; rc stay constant, indicating 
a re-isotropisation of the velocities. We tied the increase in (j) to 
the continuous radial orbit instability which always pushes a pop¬ 
ulation with low (j) / j ckc back towards the SE limit (Section 3.3 1 . 


The consequence is that - in realistically triaxial halos - the final 
dark matter core size will be chiefly determined by the total energy 
lost from baryons to dark matter, with little sensitivity to details of 
the coupling mechanism. 

Although the PG12 model can be completed neatly in this 
way, other analytic calculations or simulations based on spheri¬ 
cal symmetry will need to be considered on a case-by-case basis. 
In our exactly-spherical test cases (Figure [8]), core development is 
substantially suppressed. This certainly cautions against taking the 
results of purely spherical analyses at face value. On the other hand, 
any slight asphericity is normally sufficient to prevent the unphysi¬ 
cal angular-momentum lock-up - in particular, we verified that the 
simulations in Teyssier et al. ( [2013[ > were unaffected by this issue 
because their baryons settle into a flattened disk. The analytic re¬ 
sult is generic, so the exact shape and strength of the asphericity 
is a secondary effect in determining the final spherically-averaged 
distribution function /q. 

That said, to understand the way in which /o actually evolves 
towards the SE limit (and perhaps freezes out before it gets there) 
requires going to second order in perturbation theory, as shown in 
Appendix [B] At this point it may also become important to incor¬ 
porate self-gravity, i.e. the instantaneous connection between Sf 
and SH. This approach has been investigated more fully elsewhere, 
leading to a different set of insights regarding the onset of the radial 
orbit instability (e.g. |Saha| 199l) > as opposed to its end state. The 
present work actually suggests that the radial orbit instability can 
be cast largely as a kinematical process, and that the self-gravity 
is a secondary aspect; it would be interesting to further understand 
how these two views relate. But of more immediate practical impor¬ 
tance is to apply the broader insights about dwarf spheroidal stellar 
equilibria to observational data, something that we will attempt in 
the near future. 


ACKNOWLEDGMENTS 

All simulation analysis made use of the pynbody suite (Pontzen 
|et al.||20T3| >. AP acknowledges helpful conversations with James 
Binney, Simon White, Chervin Laporte and Chiara Tonini and sup¬ 
port from the Royal Society and, during 2013 when a significant 
part of this work was undertaken, the Oxford Martin School. JIR 
acknowledges support from SNF grant PP00P2_128540/1. FG ac¬ 
knowledges support from HST GO-1125, NSF AST-0908499. NR 
is supported by STFC and the European Research Council un¬ 
der the European Community’s Seventh Framework Programme 
(FP7/2007-2013) / ERC grant agreement no 306478-CosmicDawn. 
JD’s research is supported by Adrian Beecroft and the Oxford Mar¬ 
tin School. This research used the DiRAC Facility, jointly funded 
by STFC and the Large Facilities Capital Fund of BIS. 


References 

Adams F. C., Bloch A. M., Butler S. C., Druce J. M., Ketchum 
J. A., 2007, ApJ, 670, 1027, 0708.3101 
Arnold V. I., 1963, Russian Mathematical Surveys, 18, 9 
Barnes E. I., Lanzel P. A„ Williams L. L. R„ 2009, ApJ, 704, 372, 
0908.3873 

Barnes J., Hut P., Goodman J., 1986, ApJ, 300, 112 
Bellovary J. M., Dalcanton J. J., Babul A., Quinn T. R., Maas 
R. W., Austin C. G., Williams L. L. R., Barnes E. I., 2008, ApJ, 
685, 739, 0806.3434 


© 2015 RAS. MNRAS 000, 000-000 


















Milking the spherical cow 11 


Binney J. J., Tremaine S., 2008, Galactic dynamics. 2nd ed. 

Princeton, NJ, Princeton University Press 
Bonnivard V., Combet C., Maurin D., Walker M. G., 2015, MN- 
RAS, 446, 3002, 1407.7822 

Bryan S. E., Mao S., Kay S. T., Schaye J., Dalla Vecchia C., Booth 
C. M„ 2012, MNRAS, 422, 1863, 1109.4612 
Charbonnier A., Combet C., Daniel M.. Funk S., Hinton J. A., 
Maurin D., Power C., Read J. I., Sarkar S., Walker M. G., Wilkin¬ 
son M. I., 2011. MNRAS, 418, 1526, 1104.0412 
Cole S., Lacey C., 1996, MNRAS, 281, 716, arXiv:astro- 
ph/9510147 

Corless V. L„ King L. J., 2007, MNRAS, 380, 149, astro- 
ph/0611913 

de Zeeuw T., Merritt D., 1983, ApJ, 267, 571 
Dehnen W„ 1993, MNRAS, 265, 250 

Di Cintio A., Brook C. B., Maccio A. V., Stinson G. S., Knebe A., 
Dutton A. A., Wadsley J., 2014, MNRAS, 437, 415, 1306.0898 
Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496 
Garrison-Kimmel S., Rocha M.. Boylan-Kolchin M., Bullock J., 
Lally J., 2013, MNRAS, 433, 3539, 1301.3137 
Goldstein H.. Poole C., Safko J.. 2002, Classical Mechanics, 3rd 
edition. Addison Wesley 

Govemato F., Zolotov A., Pontzen A., Christensen C., Oh S. H., 
Brooks A. M„ Quinn T„ Shen S., Wadsley J., 2012, MNRAS, 
422, 1231, 1202.0554 

Hayashi E„ Navarro J. F„ 2006, MNRAS, 373, 1117, astro- 
ph/0608376 

Hayashi E„ Navarro J. F„ Springel V., 2007, MNRAS, 377, 50, 
arXiv: astro-ph/0612327 
Henon M„ 1973, A&A, 24, 229 
Hernquist L„ 1990, ApJ, 356, 359 

Holley-Bockelmann K., Mihos J. C., Sigurdsson S., Hernquist L., 
2001, ApJ, 549, 862, arXiv:astro-ph/0011504 
Huss A., Jain B., Steinmetz M., 1999, ApJ, 517, 64, arXiv:astro- 
ph/9803117 

Jing Y. P„ Suto Y„ 2002, ApJ, 574, 538, arXiv:astro-ph/0202064 
Kasun S. F., Evrard A. E., 2005, ApJ, 629, 781, arXiv:astro- 
ph/0408056 

Kirk D., 1992, Graphics Gems III. AP Professional, Cambridge, 
MA 

Kolmogorov A. N„ 1954, in Dokl. Akad. Nauk SSSR, Vol. 98, pp. 
527-530 

Lichtenberg A. J., Lieberman M. A., 1992, Regular and stochastic 
motion, 2nd ed 

Loebman S. R., Ivezic Z., Quinn T. R., Govemato F., Brooks 
A. M., Christensen C. R.. Juric M., 2012, ApJ. 758, L23, 
1209.2708 

Lynden-Bell D„ 1967, MNRAS. 136, 101 
MacMillan J. D., Widrow L. M., Henriksen R. N., 2006, ApJ, 653, 
43, arXiv :astro-ph/0604418 
Marechal L., Perez J., 2009, ArXiv e-prints, 0910.5177 
Martizzi D., Teyssier R., Moore B., Wentz T., 2012, MNRAS, 
422, 3081, 1112.2752 

Mashchenko S., Couchman H. M. P, Wadsley J., 2006, Nature, 
442, 539, arXiv:astro-ph/0605672 
Merritt D„ Poon M. Y„ 2004, ApJ, 606, 788, astro-ph/0302296 
Merritt D., Valluri M., 1996, ApJ, 471, 82, arXiv:astro- 
ph/9602079 

Meza A., Zamorano N., 1997, ApJ, 490, 136, astro-ph/9707004 
Moser J., 1962. Nachr. Akad. Wiss. Gttingen Math.-Phys. Kl. II, 
1 

Navarro J. F„ Frenk C. S., White S. D. M„ 1996, ApJ, 462, 563 


—, 1997, ApJ, 490, 493 

Navarro J. F., Ludlow A., Springel V., Wang J., Vogelsberger M., 
White S. D. M„ Jenkins A., Frenk C. S., Helmi A., 2010, MN¬ 
RAS, 402,21,0810.1522 

Onorbe J., Boylan-Kolchin M., Bullock J. S., Hopkins P. F., Keres 
D., Faucher-Giguere C.-A.. Quataert E., Murray N., 2015, MN¬ 
RAS, submitted, arXiv: 1502.02036 
Penarrubia J., Pontzen A., Walker M. G., Koposov S. E., 2012, 
ApJ, 759, L42, 1207.2772 

Pontzen A., Govemato F„ 2012, MNRAS, 421, 3464, 1106.0499 
—, 2013, MNRAS, 430, 121, 1210.1849 
—, 2014, Nature, 506, 171, 1402.1764 

Pontzen A., Roskar R., Stinson G., Woods R., 2013, pynbody: 
N-Body/SPH analysis for python. Astrophysics Source Code Li¬ 
brary ascl: 1305.002 

Read J. I„ Gilmore G„ 2005, MNRAS, 356, 107, arXiv:astro- 
ph/0409565 

Read J. I., Wilkinson M. I., Evans N. W., Gilmore G., Kleyna J. T., 
2006, MNRAS, 367, 387, astro-ph/0511759 
Saha P, 1991, MNRAS, 248, 494 
Saha P, Read J. I., 2009, ApJ, 690, 154, 0807.4737 
Schneider M. D., Frenk C. S., Cole S., 2012, J. Cosmology As- 
tropart. Phys., 5, 30, 1111.5616 

Stadel J.. Potter D., Moore B., Diemand J., Madau P., Zemp M., 
Kuhlen M.. Quilis V., 2009, MNRAS, 398, L21, 0808.2981 
Taylor J. E„ Navarro J. F„ 2001, ApJ, 563, 483 
Teyssier R., 2002, A&A, 385, 337, arXiv:astro-ph/0111367 
Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 
429, 3068, 1206.4895 

Tonini C., Lapi A., Salucci P, 2006, ApJ, 649, 591, arXiv:astro- 
ph/0603051 

Valluri M., Debattista V. P, Quinn T., Moore B., 2010, MNRAS, 
403, 525, 0906.4784 

van Albada T. S„ 1982. MNRAS, 201, 939 
Weinberg M. D„ 1991, ApJ, 368, 66 

Zolotov A.. Brooks A. M., Willman B., Govemato F., Pontzen A., 
Christensen C., Dekel A., Quinn T., Shen S., Wadsley J., 2012, 
ApJ, 761,71, 1207.0007 

Zorzi A. F„ Muzzio J. C., 2012, MNRAS. 423, 1955, 1204.5428 


APPENDIX A: BACKGROUND 
A1 A brief review of actions 


This Appendix contains a very brief review of actions which are 
necessary for deriving the main result in the paper. For more 
complete introductions see |Binney & Tre maine < 2008J or |GokL] 
|stein et aLl|2002| . We consider any mechanical problem described 
by phase-space coordinates q and momenta p , with Hamiltonian 
H(p, q) so that the equations of motion are 


P= ~ 


dH 

dq' 



(Al) 


where q = dq/dt and t denotes time. 

The actions can be defined starting from any coordinate sys¬ 
tem in which the motion is separately periodic in each dimension 
(i.e. for each i, qi and p ; repeat every At,). The actions Jj are then 
given by 


1 [ A <> 

J, = — / pjcij df (no sum over i). (A2) 

2 n Jo 


© 2015 RAS, MNRAS 000, 000-000 










12 A. Pontzen et al 



Figure Al. An example orbit in a 3D anisotropic harmonic oscillator. (Left) the solid line shows a portion of the orbit projected into the x-y plane. The 
equivalent orbit in an isotropic potential is closed and given by the blue dashed line. (Right) the orbit analysed in terms of J, and j, the actions for the spherical 
background Hamiltonian. Because A H is bounded (see text) the particle has to stay within a narrow band in the J r -j plane, corresponding to the linear theory 
resonance (dotted line). The blue dot shows the action for the orbit in the spherical potential (for which the actions are by definition fixed). 


By construction the actions Jj do not change under time evolution 
and are therefore integrals of motion. 

We can complete the set of 6 phase-space coordinates with 
angles 0 in such a way that equations of motion of the form {AT} 
apply. Since Jj is constant, we must then have 

ffr < A3 » 

where the first equation establishes that H can have no 0 depen¬ 
dence, and the second that consequently the 0/ each increase at a 
constant rate in time specified by the frequency £2;( J). The conve¬ 
nience of this set of coordinates is that all the time evolution of a 
particle trajectory is represented in a very simple way: 

Jj = constant, 0/ = constant + £2,f. (A4) 


Furthermore the equations of motion are canonical, which is suf¬ 
ficient to demonstrate that the coordinates are canonical - in other 
words, the measure appearing in phase space integrals is dVd 3 ©. 

We now specialise to the spherical case, with polar coordinates 
q = (r,6,<j>) and momenta p = ( f , r 2 0,r 2 sin 2 6<j>). Consider first 

1 r&t , , \ r lK -to 

Ja, = — / dttjrr" sin“ 0 = — / dtjxpr" sin - 0 = j z (A5) 

2 7i Jo 2 k Jo 

where j z is the z component of the specific angular momentum, 
j z = <j)r 2 sin 2 0 which is a constant of motion. Next consider Jg, 


Je 


1 M/e 

2k Jo 


d t0 2 r 2 




(A6) 


where j 2 = j 2 + j 2 + j 2 is the square of the total specific angular 
momentum, and 0 varies between 0 fl and 8g over the course of 
an orbit. To evaluate the integral requires a change of variables to 


relate d a and Qg to the inclination of the orbit and so to j and j z 
(e.g. [Goldstein et al.|2002| eq 10.135); the final result is that 

Je=j~iz■ (A7) 

Because linear combinations of actions are still actions ( i.e. they 
still satisfy equation ( |A3| >) one can take j z as the first and j as the 
second action in place of J<f> and Jg. 

The last action is the radial action, 

J r= 2 Io dt ' 2 = \ [. dr \/ E ~ j 2 / 2 '' 2 -®( r ). (A8) 

where r has been evaluated by energy conservation, and r librates 
between r a and rg over the period of an orbit. 

Although the complexity of expression ( | A8| > appears to make 
using the actions cumbersome, the great simplification it brings 
to the equations of motion in the background (i.e. equation [A3} 
makes the perturbation theory tractable. For that reason we have 
used the action-angle coordinates in our analytic derivation, Sec- 
tion[2] but avoided them when discussing results from simulations 
in Section[3] 


A2 Perturbed trajectories in the harmonic oscillator case 

Section [2] used Hamiltonian perturbation theory to discuss the be¬ 
haviour of particles in aspherical potentials. To connect this more 
firmly with the equations of motion, it can be helpful to study orbits 
in a specific potential and connect the solutions with the more gen¬ 
eral statements made by the perturbation theory. In this Appendix, 
we use the harmonic oscillator as such an illustrative example. 

Consider the Hamiltonian for a single particle in an 
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anisotropic harmonic oscillator potential, 

H = j + n y + n z ) + (-*“ + (1 + £ )y " + (1 — ^)z 2 ) • 

(A9) 

where K x , Tty, 7t z denote the momentum in the three cartesian direc¬ 
tions and without loss of generality e ^ 0 and <5^0. 

We want to compare the motion in this potential to the be¬ 
haviour in a sphericalised version with Hamiltonian Hq\ the change 
in the Hamiltonian between the true and sphericalised cases is given 
by 


8H = H-Hq= 1 -(oI [ey 2 - Sr) . (A10) 

We can immediately read off the first result, which is that the mag¬ 
nitude of SH is bounded. The true solution moves around on a fixed 
H surface, meaning the fractional error takes a maximum value 
given by 


M 

H 


< max(e, S). 


(All) 


This is the equivalent of the general statement that Hq variations 
are small, given by equation 

On the other hand, the angular momentum does change sig¬ 
nificantly. We can see this as follows: each of x, y and z undergoes 
oscillation at the frequencies too, (1+ £) 1//2 fi>o an d (1 — re¬ 

spectively. Assuming these new frequencies are not commensurate, 
the relative phases between the different oscillations slowly shifts 
until at some point all three separated oscillators reach x = 0, y = 0 
and z = 0 at the same moment. At this point, since the velocity re¬ 
mains finite, the angular momentum has become zero. It may take a 
number of dynamical times before this happens, but (for example) 
at the centre of dark matter halos the dynamical time is very short 
compared to the Hubble time so angular momentum conservation 
is effectively destroyed. 

All the above is illustrated in Figure |AT] The left panel con¬ 
trasts the orbits for a spherical harmonic oscillator (dashed line) and 
an aspherical oscillator (solid line) projected in the (x,y) plane. The 
latter obeys equation jA9} with e = 8 = 0.1 (and the former with 
e = 8 = 0). The orbits for the spherical case are closed because 
the frequencies are identical, so the relative phase of the x and y 
part of the motion remains fixed. Orbits in the aspherical potential 
are more complex as the relative phase of the cartesian components 
gradually changes; in fact, a particle will sometimes plunge through 
the centre of the potential. This is known as a ‘box orbit'. In more 
general triaxial potentials, a variety of orbit types are possible (e.g. 
|Merritt & Valluri| 1996| ; the importance of these will be considered 
momentarily. 

The same portion of the orbit is illustrated in the right panel, 
but now projected into the spherical actions plane (, J r ,j ). For the 
spherical case, J r and j are exactly conserved by construction and 
the orbit appears as a single point. In the aspherical case (solid line), 
J r and j are not even approximately conserved over more than a 
dynamical time. However, even then, the orbit remains close to the 
dashed line of constant Ho, as required by equation {ATT} and more 
generally by equation j8|. From this figure, it is intuitively plausible 
that the particle is equally likely to be found anywhere along the 
constant Hq contour, which is the essential result of Section [T2| 

Another way to look at the effect is to consider the relation¬ 
ship between the angular momentum and the harmonic oscillator’s 
cartesian actions, J x , J y and J z which remain constant even in the 
aspherical case. Specialising for simplicity to the case that J z = 0, 



Figure A2. Orbits in the action space of the equilibrium ‘Dwarf’ cosmo¬ 
logical halo. Dotted contours of contant spherical specific energy Hq are 
spaced at 600km 2 s -2 . (For this halo ~ 3 100km 2 s~ 2 .) The straight¬ 
ness of these contours is part of the reason that orbits can efficiently explore 
the space, as described in the text. 


one can show that 


j = 2 y/JJy sin)©* - 0 V ), (A 12) 

where 0* and 0 V are the angles conjugate to J x and J y . Only if 
0* — 0 V remains constant is j an integral of motion; as soon as the 
oscillator is aspherical, one has 

j = 2-y/Tr/y sin(0o + (£l x - £2y)0 (A13) 

so that j oscillates on the timescale 2n{£l x — £2 V ) 1 — 7T£ 1 ■ The 
situation in 3D is qualitatively similar. 

The harmonic oscillator orbits discussed above are all box or¬ 
bits. More general triaxial potentials support other orbit types (loop 
orbits, for example, which are much more tightly constrained; or 
chaotic orbits, which are even less tightly constrained than box or¬ 
bits). However the Hamiltonian analysis in the main text is general 
for all these types of possible orbit. The fraction of different orbit 
types will determine how fast and how far particles diffuse along 
the contour. In realistic cosmological dark matter halos, most or¬ 
bits in the central regions are indeed of the box or chaotic type 
i jZorzi & Muzzio|2012| l. Even with the baryonic contribution to the 
gravitational force included (which partially sphericalises the cen¬ 
tral potential), the large majority of particles remain on the same 
class of orbit as the dark progenitor jValluri et al.||2010| so long 
as feedback is strong enough to prevent long-lived central baryon 
concentrations developing (|Bryan et al.|2012). 


A3 The action space of dark matter halos 

The action-angle space of our equilibrium ‘Dwarf’ dark-matter- 
only halo is illustrated in Figure |A2] along with some particle orbits 
(solid lines) over 1.5 Gyr ~ 3fdyn- The Hamiltonian is a function 
of J r and j only for any spherical potential, and so we have sup¬ 
pressed the third action j z . As expected, the particles explore the 
space, approximately running along the Ho contours, giving rise to 
the continuous radial orbit instability described in Section |3~3l The 
contours of Hq give a great deal of dynamical information, because 
the background frequencies are defined by = dHo/dJ. These 
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frequencies are at the core of perturbation theory through equation 
& but they also determine the higher-order behaviour as follows. 

The first step to understanding the behaviour of resonances is 
to move to secular perturbation theory (e.g. section 2.4 of |Licht- | 
enberg & Lieberman 1992). Secular analysis splits resonant or¬ 
bits into two classes, known as accidentally- and intrinsically- 
degenerate. The intrinsic case refers to the situation where the res¬ 
onance condition applies globally; in other words, that n • fio is 
near-constant along lines of constant Hg. Suppose, conversely, that 
S7o did vary along these directions; then, since fig is defined by the 
normal to the Hg contours (T2g = dHg/dJ), this would imply a sig¬ 
nificant curvature of the dotted lines in Figure [A2| Since there is no 
such curvature, we can read off that the frequencies do not change 
and the dynamics is in the intrinsically-degenerate regime, giving 
rise to large-scale migrations. (The same property also means that 
the approximation £l r = £l r (Hg) used in reaching equation fH 
will be extremely accurate.) We have established using numerical 
investigations that this intrinsic degeneracy property is generic to 
any spherical action space with smooth potentials, rather than be¬ 
ing specific to cosmological halos. 


APPENDIX B: THE EVOLUTION OF F 0 

In the main text (Section |2j, we showed that a distribution func¬ 
tion / is maximally stable to linear perturbations if its sphericalised 
part /o appears ergodic, fg = fg(Hg). However we did not discuss 
the actual evolution of fg to see whether this limit is likely to be 
achieved. 

This requires time-dependent perturbation theory. We start by 
writing the collisionless Boltzmann equation, 

§£ = [HJ\ = [Hg,fo] + [Ho,8f} + [8HJo} + [8H,8f], (Bl) 

which is an exact expression. The first term vanishes identically; the 
second and third terms are linear order, and the final term is second 
order. The evolution of fg is given by taking the time derivative of 
equation and interchanging the derivative and integral: 

f-gkA* 3 ® '"•'1 

= p^Ej- j0 K^ e .«/-'- 9 ]. <B2) 

where the two linear-order terms have vanished after integrating 
over 0. Expanding the Poisson bracket and integrating the remain¬ 
ing term gives 

^=-i'En-jj(8f n 8H- n ). (B3) 

This shows that the evolution of fg is a fundamentally non-linear 
phenomenon. To make further progress, we can eliminate 8f n , 
showing that fg evolution depends only on 8H n . First, Fourier 
transform the time-dependence of fg and Sf, so that 

fg{J,t) = J dcoe ia> ' fg(J,co) and (B4) 

8f(J,B,t) = £ [dco8f n (J , co)e im+i& ' n . (B5) 

n J 

For simplicity, we will restrict attention to the case where 8H is 
constant in time. Then the evolution of 8f is given by the linear- 
order part of equation m Fourier transformed to give: 

(0) + Tig • n) 8f n = \jfj n 


) 8H n . (B6) 


This is just the time-dependent, Fourier-transformed version of 
equation G3- Note that, for 8H and 8f to be real, the Fourier 
coefficients must satisfy 

SfnM* = 8f- n (-a>) and 8H* n = 8H- n . (B7) 


These requirements are consistent with the relation |B6j. 

We can now substitute the linear-order solution fB6] for 8fn 
in equation ( |B3| > to get the leading-order evolution equation for fg: 

8H n 8H_ n n -dfg/d J1 

-—- . (do) 

CO + SIq • n 


COfg = -£n- 


dj 


Using the reality condition, equation we can pair up negative 
and positive n modes, giving an alternative version of the expres¬ 
sion that is more explicitly symmetric: 


/o = I 

n 



\SHn\ 2 n- dfg/d J 
|f2o • n| 2 — co 2 


(B9) 


where we have divided both sides by co and so the result is tech¬ 
nically only applicable for Provided the evolution of fg is 

smooth, its Fourier transform fg is continuous, so this is not a prob¬ 
lem in practice. 

Equation ( |B9| > is enough to give some insight into the relax¬ 
ation process. For simplicity, consider a perturbation consisting of 
a single, resonant ra^-mode 8H nL with n±_ ■ fig = 0. Then we can 
explicitly invert the Fourier transform to yield 

=^-jj[\8H n A 2 ri L -dfg/d J ]. (BIO) 


This is a wave equation in the nj_ direction with varying speed- 
of-sound proportional to \SH n± \. Any unevenness in the resonant 
directions will flow away at a speed proportional to the strength 
of the asphericity, showing explicitly that fg evolves towards the 
spherically ergodic limit. We hope to investigate the full behaviour 
of equation fB9} for a variety of regimes in future work. 
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